home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 050 / disk_425.arc / SIMP1.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-04-03  |  1.4 KB  |  64 lines
  1. program simp1;        { -> 273 }
  2. { integration by Simpson's method }
  3.  
  4. const    tol        = 1.0E-4;
  5. var    sum,upper,lower    : real;
  6.  
  7. external procedure cls;
  8.  
  9. function fx(x: real): real;
  10. { find f(x)=1/x }
  11. { watch out for x=0 }
  12.  
  13. begin
  14.   fx:=1.0/x
  15. end;    { function fx }
  16.  
  17. procedure simps(function f(x: real): real;
  18.         lower,upper,tol    : real;
  19.         var sum        : real);
  20.  
  21. { numerical integration by Simpson's rule }
  22. { function is f (as paramater), limits are lower and upper }
  23. { with number of regions equal to pieces }
  24. { partition is delta_x, answer is sum }
  25.  
  26. var    i        : integer;
  27.     x,delta_x,even_sum,
  28.     odd_sum,end_sum,
  29.     sum1        : real;
  30.     pieces        : integer;
  31. begin
  32.   pieces:=2;
  33.   delta_x:=(upper-lower)/pieces;
  34.   odd_sum:=f(lower+delta_x);
  35.   even_sum:=0.0;
  36.   end_sum:=f(lower)+f(upper);
  37.   sum:=(end_sum+4.0*odd_sum)*delta_x/3.0;
  38.   writeln(pieces:5,sum);
  39.   repeat
  40.     pieces:=pieces*2;
  41.     sum1:=sum;
  42.     delta_x:=(upper-lower)/pieces;
  43.     even_sum:=even_sum+odd_sum;
  44.     odd_sum:=0.0;
  45.     for i:=1 to pieces div 2 do
  46.       begin
  47.     x:=lower+delta_x*(2.0*i-1.0);
  48.     odd_sum:=odd_sum+f(x)
  49.       end;
  50.     sum:=(end_sum+4.0*odd_sum+2.0*even_sum)*delta_x/3.0;
  51.     writeln(pieces:5,sum)
  52.   until abs(sum-sum1)<=abs(tol*sum)
  53. end;    { simps }
  54.  
  55. begin        { main program }
  56.   cls;
  57.   lower:=1.0;
  58.   upper:=9.0;
  59.   writeln;
  60.   simps(fx,lower,upper,tol,sum);
  61.   writeln;
  62.   writeln(chr(7),'area= ',sum)
  63. end.
  64.